Diversity_statistics

Purpose
Required packages and data
Functions
Genotypic Diversity
- Confidence intervals for genotypic diversity
Allelic and Genotypic diversity per locus
- Tables
AMOVA
Pairwise Heatmaps
- By region
- By year

Purpose

This document will analyze and display the stastistics of diversity at the locus and genotype level.

Required packages and data

library(PramCurry)

## Loading required package: poppr
## Loading required package: adegenet
## Loading required package: ade4
##    ==========================
##     adegenet 1.4-2 is loaded
##    ==========================
## 
##  - to start, type '?adegenet'
##  - to browse adegenet website, type 'adegenetWeb()'
##  - to post questions/comments: adegenet-forum@lists.r-forge.r-project.org
## 
## 
## This is poppr version 1.1.2.99.70. To get started, type package?poppr

library(reshape2)
library(ggplot2)
library(poppr)
library(adegenet)
library(mmod)

## Loading required package: pegas
## Loading required package: ape
## 
## Attaching package: 'pegas'
## 
## The following object is masked from 'package:ape':
## 
##     mst
## 
## The following object is masked from 'package:ade4':
## 
##     amova

library(dplyr)

## 
## Attaching package: 'dplyr'
## 
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## 
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(tidyr)
options(stringsAsFactors = FALSE)
data(ramdat)
data(pop_data)
data(myPal)
sessionInfo()

## R version 3.1.2 (2014-10-31)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] tidyr_0.1         dplyr_0.2         mmod_1.2.1       
##  [4] pegas_0.6         ape_3.1-4.5       ggplot2_1.0.0    
##  [7] reshape2_1.4      PramCurry_1.0     poppr_1.1.2.99-70
## [10] adegenet_1.4-2    ade4_1.6-2       
## 
## loaded via a namespace (and not attached):
##  [1] assertthat_0.1      bitops_1.0-6        boot_1.3-11        
##  [4] caTools_1.17.1      colorspace_1.2-4    digest_0.6.4       
##  [7] evaluate_0.5.5      fastmatch_1.0-4     formatR_1.0        
## [10] gdata_2.13.3        ggmap_2.3           grid_3.1.2         
## [13] gtable_0.1.2        gtools_3.4.1        htmltools_0.2.6    
## [16] httpuv_1.3.0        igraph_0.7.1        knitr_1.6          
## [19] labtools_0.4.1      lattice_0.20-29     lubridate_1.3.3    
## [22] mapproj_1.2-2       maps_2.3-9          MASS_7.3-34        
## [25] Matrix_1.1-4        memoise_0.2.1       munsell_0.4.2      
## [28] nlme_3.1-117        parallel_3.1.2      permute_0.8-3      
## [31] phangorn_1.99-7     plotrix_3.5-7       plyr_1.8.1         
## [34] png_0.1-7           proto_0.3-10        Rcpp_0.11.2        
## [37] RgoogleMaps_1.2.0.6 rjson_0.2.14        RJSONIO_1.3-0      
## [40] rmarkdown_0.3.10    scales_0.2.4        shiny_0.10.1       
## [43] stringr_0.6.2       tools_3.1.2         vegan_2.0-10       
## [46] xtable_1.7-4        yaml_2.1.13

myTheme <- theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))

Functions

get_stats <- function(z){
  mat <- matrix(nrow = nrow(z), ncol = 4, 
                dimnames = list(Pop = rownames(z), 
                                Index = c("H", "G", "Hexp", "E.5")
                                )
                )
  N     <- rowSums(z)
  H     <- vegan::diversity(z)
  G     <- vegan::diversity(z, "inv")
  Simp  <- vegan::diversity(z, "simp")
  nei   <- (N/(N-1)) * Simp
  E.5   <- (G - 1)/(exp(H) -1)
  mat[] <- c(H, G, nei, E.5)
  return(mat)
}

boot_stats <- function(x, i){
  res        <- numeric(4)
  names(res) <- c("H", "G", "Hexp", "E.5")
  z     <- table(x[i])
  N     <- sum(z)
  H     <- vegan::diversity(z)
  G     <- vegan::diversity(z, "inv")
  Simp  <- vegan::diversity(z, "simp")
  nei   <- (N/(N-1)) * Simp
  E.5   <- (G - 1)/(exp(H) -1)
  res[] <- c(H, G, nei, E.5)
  return(res)
}

extract_samples <- function(x) rep(1:length(x), x)

do_boot <- function(tab, n, ...){
  res <- apply(tab, 1, function(x) boot::boot(extract_samples(x), boot_stats, n, ...))
  return(res)
}

get_ci <- function(x, lb, ub){
  res <- apply(x$t, 2, quantile, c(lb, ub), na.rm = TRUE)
  return(res)
}

get_all_ci <- function(res, ci = 95){
  lower_bound  <- (100 - ci)/200
  upper_bound  <- 1 - lower_bound
  funval       <- matrix(numeric(8), nrow = 2)
  CI           <- vapply(res, FUN = get_ci, FUN.VALUE = funval, 
                         lower_bound, upper_bound)
  dCI          <- dimnames(CI)
  dimnames(CI) <- list(CI    = dCI[[1]], 
                       Index = c("H", "G", "Hexp", "E.5"),
                       Pop   = dCI[[3]])
  
  return(CI)
}


boot_ci <- function(tab, n = 1000, ci = 95, total = TRUE, ...){
  if (!is.matrix(tab) & is.genind(tab)){
    tab <- mlg.table(tab, total = total, bar = FALSE)
  }
  res  <- do_boot(tab, n, ...)
  orig <- get_stats(tab)
  orig <- melt(orig)
  CI   <- get_all_ci(res, ci = ci)
  samp <- vapply(res, "[[", FUN.VALUE = res[[1]]$t, "t")
  dimnames(samp) <- list(NULL, 
                         Index = c("H", "G", "Hexp", "E.5"),
                         Pop = rownames(tab))
  pl <- ggplot(melt(samp), aes(x = factor(Pop), y = value, group = Pop)) + 
    geom_boxplot() + 
    geom_point(aes(color = factor(Pop), x = factor(Pop), y = value), 
               size = 5, pch = 16, data = orig) +
    xlab("Population") + labs(color = "Observed") +
    facet_wrap(~Index, scales = "free_y") + myTheme
  print(pl)
  return(CI)
}

calc_loc_table <- function(dat, hier, combine = FALSE, lev = "allele", cc = TRUE){
  if (cc){
    dat <- clonecorrect(dat, hier = hier, combine = combine)
  }
  theTable <- locus_table(dat, information = FALSE, lev = lev)
  theHier <- gethierarchy(dat, hier, combine = combine)
  if (combine){
    setpop(dat) <- hier
    theHier <- levels(theHier[[length(theHier)]])
  } else {
    setpop(dat) <- as.formula(paste0("~", names(theHier)[1]))
    theHier <- levels(theHier[[1]])
  }
  loc_tab <- lapply(theHier, 
                   function(x){
                     locus_table(dat, 
                                 population = x, 
                                 information = FALSE,
                                 lev = lev)
                   })
  names(loc_tab) <- theHier
  return(list(pops = loc_tab, total = theTable))
}


plot_loc_table <- function(tab){
#   tab <- theList[["pops"]]
#   tot <- melt(theList["total"])
  tab <- melt(tab)
  ggplot(tab, aes(y = value, x = L2, fill = L2, linetype = NA)) +
    geom_bar(stat = "identity") +
    geom_hline(aes(yintercept = value, linetype = L1), 
               data = tab[tab$L1 == "total", ]) + 
    facet_grid(summary ~ locus, scales = "free_y") +
    scale_linetype_discrete(labels = "", name = "Pooled") +
    myTheme
}



amova_pair <- function(pops, dat, hier, ...){
  dat <- popsub(dat, pops)
  if (any(table(pop(dat)) < 3)){
    return(NULL)
  }
  poppr.amova(dat, hier, ...)
}

pairwise_amova <- function(x, hier, ...){
  pops         <- x@pop.names
  pop_combs    <- combn(pops, 2)
  xlist        <- apply(pop_combs, 2, amova_pair, x, hier, ...)
  names(xlist) <- apply(pop_combs, 2, paste, collapse = " : ")
  return(xlist)
}

pairwise_amova_test <- function(pairlist, nrepet = 99){
  res <- lapply(names(pairlist), print_amova_pairs, pairlist, nrepet)
  names(res) <- names(pairlist)
  return(res)
}

print_amova_pairs <- function(pairname, pairlist, nrepet){
  thePair <- pairlist[[pairname]]
  if (!is.null(thePair)){
    cat(pairname, "\n")
    theTest <- randtest(thePair, nrepet = nrepet)
    try(plot(theTest))
    return(theTest)
  } else {
    return(NULL)
  }
}


getpval <- function(x){
  if (is.null(x$pvalue)){
    return(as.numeric(rep(NA, 4)))
  } else {
    pvals <- x$pvalue
  }
  pvals[x$rep == 0] <- NA
  return(pvals)
}

getphi <- function(x){
  if (is.null(x)) return(numeric(4))
  return(x$statphi$Phi)
}

refactor <- function(strings, factors){
  u_factors <- unique(factors)
  u_factors[match(strings, u_factors)]
}

make_amova_table <- function(am, amt){
  tot <- nrow(am$results)
  res <- data.frame(list(am$results[-tot, c("Df", "Sum Sq")], 
                         Percent = am$componentsofcovariance[-tot, 2],
                         Pval    = rev(amt$pvalue), 
                         Phi     = am$statphi$Phi[-tot]))
  res <- as.matrix(res)
  colnames(res) <- c("d.f.", "Sum of Squares", "Percent variation", "P", 
                     "Phi statistic")
  names(dimnames(res)) <- c("levels", "statistic")
  return(res)
}

make_amova_printable <- function(amtab, amtabcc){
  am_array <- array(dim      = c(dim(amtab), 2),
                    dimnames = c(dimnames(amtab), 
                                 list(c("full", "clone-corrected"))))
  am_array[, , 1] <- amtab
  am_array[, , 2] <- amtabcc
  tabfun <- function(x){
    x <- paste0(paste0(signif(x, 3), collapse = " ("), ")")
    return(x)
  }
  res <- apply(am_array, c(1, 2), tabfun)
  return(res)
}

ci2string <- function(x, sig = 3, colString = "-"){
  if (all(is.na(x))) return("(-)")
  paste0("(", paste(signif(x, sig), collapse = colString),")")
}

add_ci <- function(dtable, ciarray, sig = 3, colString = "-"){
  
  dtable <- data.frame(lapply(dtable, function(x, sig){
    if (is.numeric(x)) x <- signif(x, sig)
    return(x)
  }, sig))
  for (i in colnames(ciarray)){
    to_add <- apply(ciarray[, i, ], 2, ci2string, sig, colString)
    dtable[[i]] <- paste(signif(dtable[[i]], sig), to_add)
  }
  return(dtable)
}

Genotypic Diversity

(GDxYear      <- poppr(ramdat, 
                      sample = 999, 
                      quiet = TRUE, 
                      hist = FALSE))

##      Pop   N MLG eMLG    SE     H     G  Hexp   E.5      Ia  p.Ia    rbarD
## 1   2001  39  10 3.65 1.161 1.187  1.90 0.486 0.395 -0.0864 0.869 -0.06765
## 2   2002  36   9 4.20 1.133 1.373  2.34 0.589 0.454  0.4555 0.001  0.19332
## 3   2003 102  20 4.75 1.301 1.811  2.92 0.665 0.376  0.4262 0.001  0.14326
## 4   2004  35   8 4.56 1.007 1.575  3.53 0.738 0.661  0.5282 0.001  0.18405
## 5   2005   2   2 2.00 0.000 0.693  2.00 1.000 1.000      NA    NA       NA
## 6   2006   1   1 1.00 0.000 0.000  1.00   NaN   NaN      NA    NA       NA
## 7   2010   1   1 1.00 0.000 0.000  1.00   NaN   NaN      NA    NA       NA
## 8   2011  68   6 2.15 0.827 0.563  1.32 0.245 0.420  0.1720 0.001  0.18094
## 9   2012  20   7 5.01 0.916 1.623  4.00 0.789 0.737 -0.1417 0.869 -0.17493
## 10  2013 179  47 7.74 1.204 3.148 13.77 0.933 0.573  0.0623 0.015  0.02925
## 11  2014  30  17 8.09 1.029 2.671 12.16 0.949 0.830 -0.0112 0.503 -0.00617
## 12 Total 513  70 6.95 1.332 2.983  8.64 0.886 0.408  0.2620 0.001  0.08159
##     p.rD   File
## 1  0.999 ramdat
## 2  0.001 ramdat
## 3  0.001 ramdat
## 4  0.001 ramdat
## 5     NA ramdat
## 6     NA ramdat
## 7     NA ramdat
## 8  0.001 ramdat
## 9  1.000 ramdat
## 10 0.001 ramdat
## 11 0.555 ramdat
## 12 0.001 ramdat

(GDxRegion    <- poppr(setpop(ramdat, ~ZONE2), 
                      sample = 999, 
                      quiet = TRUE, 
                      hist = FALSE))

##          Pop   N MLG eMLG    SE     H    G  Hexp   E.5      Ia  p.Ia
## 1    JHallCr 244  30 4.89 1.335 1.972 3.13 0.683 0.345  0.3322 0.001
## 2 NFChetHigh 114  35 6.81 1.342 2.735 7.07 0.866 0.421  0.1426 0.001
## 3      Coast  34  12 5.91 1.169 2.048 5.56 0.845 0.675  0.2514 0.011
## 4   HunterCr  66   4 1.88 0.690 0.423 1.24 0.198 0.460 -0.0440 0.770
## 5   Winchuck  35   9 4.29 1.072 1.474 2.88 0.672 0.559 -0.0196 0.593
## 6 ChetcoMain  16   7 5.00 0.926 1.450 2.84 0.692 0.565  0.3748 0.017
## 7  PistolRSF   4   2 2.00 0.000 0.562 1.60 0.500 0.795  0.0000 0.459
## 8      Total 513  70 6.95 1.332 2.983 8.64 0.886 0.408  0.2620 0.001
##     rbarD  p.rD                   File
## 1  0.1031 0.001 setpop(ramdat, ~ZONE2)
## 2  0.0660 0.001 setpop(ramdat, ~ZONE2)
## 3  0.2828 0.001 setpop(ramdat, ~ZONE2)
## 4 -0.0471 1.000 setpop(ramdat, ~ZONE2)
## 5 -0.0138 0.756 setpop(ramdat, ~ZONE2)
## 6  0.3957 0.001 setpop(ramdat, ~ZONE2)
## 7     NaN    NA setpop(ramdat, ~ZONE2)
## 8  0.0816 0.001 setpop(ramdat, ~ZONE2)

(GDxYearccYR  <- poppr(ramdat, 
                      clonecorrect = TRUE, 
                      hier = ~Pop/ZONE2, 
                      sample = 999, 
                      quiet = TRUE, 
                      hist = FALSE))

##      Pop   N MLG  eMLG    SE     H    G  Hexp   E.5      Ia  p.Ia   rbarD
## 1   2001  10  10 10.00 0.000 2.303 10.0 1.000 1.000 -0.3622 0.961 -0.1891
## 2   2002   9   9  9.00 0.000 2.197  9.0 1.000 1.000  0.1957 0.207  0.0668
## 3   2003  21  20  9.79 0.410 2.979 19.2 0.995 0.974  0.2319 0.029  0.0615
## 4   2004   8   8  8.00 0.000 2.079  8.0 1.000 1.000  0.1611 0.262  0.0544
## 5   2005   2   2  2.00 0.000 0.693  2.0 1.000 1.000      NA    NA      NA
## 6   2006   1   1  1.00 0.000 0.000  1.0   NaN   NaN      NA    NA      NA
## 7   2010   1   1  1.00 0.000 0.000  1.0   NaN   NaN      NA    NA      NA
## 8   2011   6   6  6.00 0.000 1.792  6.0 1.000 1.000 -0.3023 0.783 -0.3226
## 9   2012   7   7  7.00 0.000 1.946  7.0 1.000 1.000 -0.3223 0.867 -0.3334
## 10  2013  68  47  9.42 0.710 3.724 35.6 0.986 0.855 -0.1275 0.998 -0.0508
## 11  2014  19  17  9.47 0.598 2.799 15.7 0.988 0.953 -0.1058 0.818 -0.0545
## 12 Total 152  70  9.19 0.829 3.939 38.4 0.980 0.742  0.0585 0.040  0.0155
##     p.rD   File
## 1  1.000 ramdat
## 2  0.145 ramdat
## 3  0.024 ramdat
## 4  0.212 ramdat
## 5     NA ramdat
## 6     NA ramdat
## 7     NA ramdat
## 8  1.000 ramdat
## 9  1.000 ramdat
## 10 1.000 ramdat
## 11 0.978 ramdat
## 12 0.035 ramdat

(GDxRegionccR <- poppr(ramdat, 
                      clonecorrect = TRUE, 
                      hier = ~ZONE2, 
                      sample = 999, 
                      quiet = TRUE, 
                      hist = FALSE))

##          Pop  N MLG eMLG       SE     H    G  Hexp   E.5      Ia  p.Ia
## 1    JHallCr 30  30 10.0      NaN 3.401 30.0 1.000 1.000  0.0537 0.251
## 2 NFChetHigh 35  35 10.0 8.89e-07 3.555 35.0 1.000 1.000 -0.2156 1.000
## 3      Coast 12  12 10.0      NaN 2.485 12.0 1.000 1.000  0.0248 0.388
## 4   HunterCr  4   4  4.0 0.00e+00 1.386  4.0 1.000 1.000 -0.5714 0.855
## 5   Winchuck  9   9  9.0 0.00e+00 2.197  9.0 1.000 1.000 -0.3824 0.961
## 6 ChetcoMain  7   7  7.0 0.00e+00 1.946  7.0 1.000 1.000 -0.0945 0.581
## 7  PistolRSF  2   2  2.0 0.00e+00 0.693  2.0 1.000 1.000      NA    NA
## 8      Total 99  70  9.6 6.06e-01 4.113 51.9 0.991 0.846  0.0409 0.146
##     rbarD  p.rD   File
## 1  0.0136 0.250 ramdat
## 2 -0.0803 1.000 ramdat
## 3  0.0249 0.239 ramdat
## 4 -0.5774 0.999 ramdat
## 5 -0.1954 1.000 ramdat
## 6 -0.0988 0.924 ramdat
## 7      NA    NA ramdat
## 8  0.0107 0.135 ramdat

(Totcc        <- poppr(ramdat, 
                      clonecorrect = TRUE, 
                      hier = NA, 
                      sublist = "Total",
                      sample = 999, 
                      quiet = TRUE, 
                      hist = FALSE))

##     Pop  N MLG eMLG SE    H  G Hexp E.5      Ia  p.Ia    rbarD  p.rD
## 1 Total 70  70   70  0 4.25 70    1   1 -0.0216 0.659 -0.00551 0.663
##     File
## 1 ramdat

Confidence intervals for genotypic diversity

# By Year:
set.seed(9001)
system.time(year_simp <- boot_ci(ramdat, n = 9999, total = TRUE, 
                      parallel = "multicore", ncpus = 4L))

## Warning: Removed 19998 rows containing non-finite values (stat_boxplot).

## Warning: Removed 25017 rows containing non-finite values (stat_boxplot).

## Warning: Removed 2 rows containing missing values (geom_point).

## Warning: Removed 2 rows containing missing values (geom_point).

##    user  system elapsed 
##  82.880   7.476  36.146

# By Population:
set.seed(9001)
system.time(zone_simp <- boot_ci(setpop(ramdat, ~ZONE2), n = 9999, total = TRUE, 
                      parallel = "multicore", ncpus = 4L))

## Warning: Removed 3196 rows containing non-finite values (stat_boxplot).

##    user  system elapsed 
##  61.808   6.224  26.320

write.table(add_ci(GDxYear, year_simp, colString = "-")[-nrow(GDxYear), -length(GDxYear)], 
            file = "GD.csv", 
            na = "-",
            row.names = FALSE,
            sep = ",", col.names = TRUE)
write.table(add_ci(GDxRegion, zone_simp, colString = "-")[-length(GDxYear)],
            file = "GD.csv", 
            na = "-",
            row.names = FALSE,
            sep = ",", col.names = FALSE, append = TRUE)

Allelic and Genotypic diversity per locus

(loc_whole <- locus_table(ramdat)) # Whole data set

## 
##  allele = Number of observed alleles
##  1-D = Simpson index
##  Hexp = Nei's 1978 expected heterozygosity
## ------------------------------------------

##           summary
## locus      allele   1-D  Hexp Evenness
##   PrMS6A1    2.00  0.50  1.00     1.00
##   Pr9C3A1    2.00  0.50  1.00     1.00
##   PrMS39A1   5.00  0.60  0.75     0.84
##   PrMS45A1   3.00  0.50  0.75     0.99
##   PrMS43A1  18.00  0.81  0.86     0.67
##   mean       6.00  0.58  0.87     0.90

(loc_geno  <- locus_table(clonecorrect(ramdat, hier = NA))) # By genotype

## 
##  allele = Number of observed alleles
##  1-D = Simpson index
##  Hexp = Nei's 1978 expected heterozygosity
## ------------------------------------------

##           summary
## locus      allele   1-D  Hexp Evenness
##   PrMS6A1    2.00  0.50  0.99     0.99
##   Pr9C3A1    2.00  0.49  0.99     0.99
##   PrMS39A1   5.00  0.62  0.78     0.78
##   PrMS45A1   3.00  0.51  0.76     0.95
##   PrMS43A1  18.00  0.89  0.95     0.78
##   mean       6.00  0.60  0.89     0.90

loc_year_allele <- calc_loc_table(ramdat, ~Pop/ZONE2)
(lya <- plot_loc_table(loc_year_allele) + 
   ggtitle("Allelic Diversity") + labs(list(x = "Year", fill = "Year")))

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

loc_year_geno <- calc_loc_table(ramdat, ~Pop/ZONE2, lev = "genotype")
(lyg <- plot_loc_table(loc_year_geno) + 
   ggtitle("Genotype Diversity") + labs(list(x = "Year", fill = "Year")))

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 7 rows containing missing values (position_stack).

## Warning: Removed 8 rows containing missing values (position_stack).

## Warning: Removed 4 rows containing missing values (position_stack).

## Warning: Removed 8 rows containing missing values (position_stack).

## Warning: Removed 3 rows containing missing values (position_stack).

## Warning: Removed 3 rows containing missing values (position_stack).

## Warning: Removed 7 rows containing missing values (position_stack).

## Warning: Removed 8 rows containing missing values (position_stack).

## Warning: Removed 4 rows containing missing values (position_stack).

## Warning: Removed 8 rows containing missing values (position_stack).

## Warning: Removed 3 rows containing missing values (position_stack).

## Warning: Removed 3 rows containing missing values (position_stack).

loc_region_allele <- calc_loc_table(ramdat, ~ZONE2)
(lya <- plot_loc_table(loc_region_allele) + 
   ggtitle("Allelic Diversity") + labs(list(x = "Region", fill = "Region")))

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

loc_region_geno <- calc_loc_table(ramdat, ~ZONE2, lev = "genotype")
(lyg <- plot_loc_table(loc_region_geno) + 
   ggtitle("Genotype Diversity") + labs(list(x = "Region", fill = "Region")))

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 6 rows containing missing values (position_stack).

## Warning: Removed 7 rows containing missing values (position_stack).

## Warning: Removed 2 rows containing missing values (position_stack).

## Warning: Removed 5 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 6 rows containing missing values (position_stack).

## Warning: Removed 7 rows containing missing values (position_stack).

## Warning: Removed 2 rows containing missing values (position_stack).

## Warning: Removed 5 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

## Warning: Removed 1 rows containing missing values (position_stack).

loc_region_allele_summary <- melt(loc_region_allele) %>% 
  filter(locus == "mean") %>% 
  mutate(region = ifelse(L1 == "pops", L2, L1)) %>% 
  select(locus, summary, value, region) %>% 
  dcast(region ~ summary) %>%
  arrange(region)
loc_year_allele_summary <- melt(loc_year_allele) %>% 
  filter(locus == "mean") %>% 
  mutate(year = ifelse(L1 == "pops", L2, L1)) %>% 
  select(locus, summary, value, year) %>% 
  dcast(year ~ summary)

Tables

Mean Allelic Diversity by Region

region	allele	1-D	Hexp	Evenness
ChetcoMain	3.0	0.5837	0.9655	0.9432
Coast	3.8	0.5944	0.9733	0.9597
HunterCr	2.6	0.5562	0.9615	0.9526
JHallCr	4.4	0.5720	0.9145	0.8908
NFChetHigh	4.8	0.5984	0.8961	0.8940
PistolRSF	2.4	0.5500	1.0000	1.0000
Winchuck	3.4	0.5630	0.9387	0.9165
total	6.0	0.6000	0.8909	0.9008

Mean Allelic Diversity by Year

year	allele	1-D	Hexp	Evenness
2001	3.6	0.5640	0.9370	0.9099
2002	3.0	0.5642	0.9443	0.9306
2003	4.2	0.5646	0.9025	0.8931
2004	2.8	0.5031	0.8441	0.8548
2005	2.4	0.5000	0.9000	0.9180
2006	2.0	0.5000	1.0000	1.0000
2010	2.0	0.5000	1.0000	1.0000
2011	3.0	0.5833	0.9767	0.9638
2012	3.2	0.5755	0.9588	0.9498
2013	5.4	0.5979	0.8918	0.8987
2014	3.8	0.5909	0.9703	0.9347
total	6.0	0.5929	0.8823	0.8976

Total Allelic Diversity

	allele	1-D	Hexp	Evenness
PrMS6A1	2.0000	0.4963	0.9927	0.9927
Pr9C3A1	2.0000	0.4950	0.9900	0.9901
PrMS39A1	5.0000	0.6229	0.7786	0.7817
PrMS45A1	3.0000	0.5058	0.7587	0.9533
PrMS43A1	18.0000	0.8932	0.9457	0.7837
mean	6.0000	0.6026	0.8931	0.9003

AMOVA

Including within individual variance

Year/Region

namehierarchy(ramdat) <- ~Year/ZONE2/ZONE1

(year_z2   <- poppr.amova(ramdat, ~Year/ZONE2))

## 
##  No missing values detected.

## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
## 
## $results
##                                Df     Sum Sq   Mean Sq
## Between Year                   10   79.79508 7.9795082
## Between ZONE2 Within Year      12   29.77109 2.4809238
## Between samples Within ZONE2  490  141.01083 0.2877772
## Within samples                513 1243.50000 2.4239766
## Total                        1025 1494.07700 1.4576361
## 
## $componentsofcovariance
##                                                Sigma          %
## Variations  Between Year                  0.02997373   2.035906
## Variations  Between ZONE2 Within Year     0.08640422   5.868836
## Variations  Between samples Within ZONE2 -1.06809970 -72.548560
## Variations  Within samples                2.42397661 164.643817
## Total variations                          1.47225486 100.000000
## 
## $statphi
##                           Phi
## Phi-samples-total -0.64643817
## Phi-samples-ZONE2 -0.78775566
## Phi-ZONE2-Year     0.05990803
## Phi-Year-total     0.02035906

(year_z2cc <- poppr.amova(ramdat, ~Year/ZONE2, clonecorrect = TRUE))

## 
##  No missing values detected.

## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
## 
## $results
##                               Df     Sum Sq   Mean Sq
## Between Year                  10  10.523948 1.0523948
## Between ZONE2 Within Year     12   9.459087 0.7882573
## Between samples Within ZONE2 129  71.609070 0.5551091
## Within samples               152 359.000000 2.3618421
## Total                        303 450.592105 1.4871027
## 
## $componentsofcovariance
##                                                 Sigma           %
## Variations  Between Year                  0.009753153   0.6553138
## Variations  Between ZONE2 Within Year     0.020089101   1.3497856
## Variations  Between samples Within ZONE2 -0.903366518 -60.6971504
## Variations  Within samples                2.361842105 158.6920509
## Total variations                          1.488317840 100.0000000
## 
## $statphi
##                            Phi
## Phi-samples-total -0.586920509
## Phi-samples-ZONE2 -0.619390908
## Phi-ZONE2-Year     0.013586893
## Phi-Year-total     0.006553138

set.seed(9001)
system.time(year_z2.t <- randtest(year_z2, nrepet = 999))

##    user  system elapsed 
## 153.881   1.489 155.397

plot(year_z2.t)

set.seed(9001)
system.time(year_z2cc.t <- randtest(year_z2cc, nrepet = 999))

##    user  system elapsed 
##  17.137   0.212  17.353

plot(year_z2cc.t)

date()

## [1] "Fri Oct 24 16:21:38 2014"

Region/Year

namehierarchy(ramdat) <- ~Year/ZONE2/ZONE1

(z2_year   <- poppr.amova(ramdat, ~ZONE2/Year))

## 
##  No missing values detected.

## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
## 
## $results
##                               Df     Sum Sq    Mean Sq
## Between ZONE2                  6   98.53739 16.4228985
## Between Year Within ZONE2     16   11.02878  0.6892985
## Between samples Within Year  490  141.01083  0.2877772
## Within samples               513 1243.50000  2.4239766
## Total                       1025 1494.07700  1.4576361
## 
## $componentsofcovariance
##                                              Sigma           %
## Variations  Between ZONE2                0.1283283   8.5740970
## Variations  Between Year Within ZONE2    0.0124923   0.8346573
## Variations  Between samples Within Year -1.0680997 -71.3637662
## Variations  Within samples               2.4239766 161.9550118
## Total variations                         1.4966975 100.0000000
## 
## $statphi
##                            Phi
## Phi-samples-total -0.619550118
## Phi-samples-Year  -0.787755656
## Phi-Year-ZONE2     0.009129331
## Phi-ZONE2-total    0.085740970

(z2_yearcc <- poppr.amova(ramdat, ~ZONE2/Year, clonecorrect = TRUE))

## 
##  No missing values detected.

## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
## 
## $results
##                              Df     Sum Sq   Mean Sq
## Between ZONE2                 6  11.582493 1.9304155
## Between Year Within ZONE2    16   8.400542 0.5250339
## Between samples Within Year 129  71.609070 0.5551091
## Within samples              152 359.000000 2.3618421
## Total                       303 450.592105 1.4871027
## 
## $componentsofcovariance
##                                                Sigma           %
## Variations  Between ZONE2                0.039160694   2.6194612
## Variations  Between Year Within ZONE2   -0.002645858  -0.1769816
## Variations  Between samples Within Year -0.903366518 -60.4262411
## Variations  Within samples               2.361842105 157.9837615
## Total variations                         1.494990423 100.0000000
## 
## $statphi
##                            Phi
## Phi-samples-total -0.579837615
## Phi-samples-Year  -0.619390908
## Phi-Year-ZONE2    -0.001817423
## Phi-ZONE2-total    0.026194612

set.seed(9001)
system.time(z2_year.t <- randtest(z2_year, nrepet = 999))

##    user  system elapsed 
## 154.778   0.348 155.143

plot(z2_year.t)

set.seed(9001)
system.time(z2_yearcc.t <- randtest(z2_yearcc, nrepet = 999))

##    user  system elapsed 
##  17.025   0.012  17.034

plot(z2_yearcc.t)

date()

## [1] "Fri Oct 24 16:42:51 2014"

Excluding within individual variance

Year/Region

namehierarchy(ramdat) <- ~Year/ZONE2/ZONE1

(WF_year_z2   <- poppr.amova(ramdat, ~Year/ZONE2, within = FALSE))

## 
##  No missing values detected.

## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
## 
## $results
##                              Df Sum Sq Mean Sq
## Between Year                 10  159.7 15.9691
## Between samples Within Year  12   59.5  4.9581
## Within samples              490  280.8  0.5731
## Total                       512  500.0  0.9766
## 
## $componentsofcovariance
##                                          Sigma      %
## Variations  Between Year                0.1203  11.58
## Variations  Between samples Within Year 0.3455  33.26
## Variations  Within samples              0.5731  55.16
## Total variations                        1.0389 100.00
## 
## $statphi
##                      Phi
## Phi-samples-total 0.4484
## Phi-samples-Year  0.3761
## Phi-Year-total    0.1158

(WF_year_z2cc <- poppr.amova(ramdat, ~Year/ZONE2, clonecorrect = TRUE, 
                          within = FALSE))

## 
##  No missing values detected.

## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
## 
## $results
##                              Df Sum Sq Mean Sq
## Between Year                 10  21.01   2.101
## Between samples Within Year  12  19.14   1.595
## Within samples              129 141.29   1.095
## Total                       151 181.44   1.202
## 
## $componentsofcovariance
##                                           Sigma       %
## Variations  Between Year                0.03654   3.001
## Variations  Between samples Within Year 0.08615   7.073
## Variations  Within samples              1.09524  89.926
## Total variations                        1.21794 100.000
## 
## $statphi
##                       Phi
## Phi-samples-total 0.10074
## Phi-samples-Year  0.07292
## Phi-Year-total    0.03001

set.seed(9001)
system.time(WF_year_z2.t <- randtest(WF_year_z2, nrepet = 9999))

##    user  system elapsed 
##   35.34    0.27   35.98

plot(WF_year_z2.t)

plot of chunk AMOVA_WF

set.seed(9001)
system.time(WF_year_z2cc.t <- randtest(WF_year_z2cc, nrepet = 9999))

##    user  system elapsed 
##  18.973   0.097  19.129

plot(WF_year_z2cc.t)

plot of chunk AMOVA_WF

date()

## [1] "Tue Nov  4 09:32:42 2014"

Region/Year

namehierarchy(ramdat) <- ~Year/ZONE2/ZONE1

(WF_z2_year   <- poppr.amova(ramdat, ~ZONE2/Year, within = FALSE))

## 
##  No missing values detected.

## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
## 
## $results
##                               Df Sum Sq Mean Sq
## Between ZONE2                  6 196.69 32.7822
## Between samples Within ZONE2  16  22.49  1.4059
## Within samples               490 280.82  0.5731
## Total                        512 500.00  0.9766
## 
## $componentsofcovariance
##                                            Sigma       %
## Variations  Between ZONE2                0.51127  44.999
## Variations  Between samples Within ZONE2 0.05182   4.561
## Variations  Within samples               0.57309  50.440
## Total variations                         1.13618 100.000
## 
## $statphi
##                       Phi
## Phi-samples-total 0.49560
## Phi-samples-ZONE2 0.08293
## Phi-ZONE2-total   0.44999

(WF_z2_yearcc <- poppr.amova(ramdat, ~ZONE2/Year, clonecorrect = TRUE, 
                          within = FALSE))

## 
##  No missing values detected.

## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
## 
## $results
##                               Df Sum Sq Mean Sq
## Between ZONE2                  6  22.95   3.825
## Between samples Within ZONE2  16  17.20   1.075
## Within samples               129 141.29   1.095
## Total                        151 181.44   1.202
## 
## $componentsofcovariance
##                                              Sigma        %
## Variations  Between ZONE2                 0.152442  12.2526
## Variations  Between samples Within ZONE2 -0.003527  -0.2835
## Variations  Within samples                1.095242  88.0308
## Total variations                          1.244157 100.0000
## 
## $statphi
##                         Phi
## Phi-samples-total  0.119692
## Phi-samples-ZONE2 -0.003231
## Phi-ZONE2-total    0.122526

set.seed(9001)
system.time(WF_z2_year.t <- randtest(WF_z2_year, nrepet = 9999))

##    user  system elapsed 
##  35.077   0.341  35.820

plot(WF_z2_year.t)

plot of chunk AMOVAzy_WF

set.seed(9001)
system.time(WF_z2_yearcc.t <- randtest(WF_z2_yearcc, nrepet = 9999))

##    user  system elapsed 
##  17.879   0.099  18.057

plot(WF_z2_yearcc.t)

plot of chunk AMOVAzy_WF

date()

## [1] "Tue Nov  4 09:33:37 2014"

Writing tables

z2year_table   <- make_amova_table(WF_z2_year, WF_z2_year.t)
z2yearcc_table <- make_amova_table(WF_z2_yearcc, WF_z2_yearcc.t)
(z2year_full_table <- make_amova_printable(z2year_table, z2yearcc_table))

##                               statistic
## levels                         d.f.        Sum of Squares
##   Between ZONE2                "6 (6)"     "197 (23)"    
##   Between samples Within ZONE2 "16 (16)"   "22.5 (17.2)" 
##   Within samples               "490 (129)" "281 (141)"   
##                               statistic
## levels                         Percent variation P              
##   Between ZONE2                "45 (12.3)"       "1e-04 (1e-04)"
##   Between samples Within ZONE2 "4.56 (-0.283)"   "1e-04 (0.446)"
##   Within samples               "50.4 (88)"       "1e-04 (1e-04)"
##                               statistic
## levels                         Phi statistic      
##   Between ZONE2                "0.496 (0.12)"     
##   Between samples Within ZONE2 "0.0829 (-0.00323)"
##   Within samples               "0.45 (0.123)"

write.table(z2year_full_table, file = "zone_by_year.csv", row.names = TRUE,
            col.names = NA, sep = ",")

yearz2_table   <- make_amova_table(WF_year_z2, WF_year_z2.t)
yearz2cc_table <- make_amova_table(WF_year_z2cc, WF_year_z2cc.t)
(yearz2_full_table <- make_amova_printable(yearz2_table, yearz2cc_table))

##                              statistic
## levels                        d.f.        Sum of Squares Percent variation
##   Between Year                "10 (10)"   "160 (21)"     "11.6 (3)"       
##   Between samples Within Year "12 (12)"   "59.5 (19.1)"  "33.3 (7.07)"    
##   Within samples              "490 (129)" "281 (141)"    "55.2 (89.9)"    
##                              statistic
## levels                        P               Phi statistic   
##   Between Year                "0.366 (0.175)" "0.448 (0.101)" 
##   Between samples Within Year "1e-04 (2e-04)" "0.376 (0.0729)"
##   Within samples              "1e-04 (1e-04)" "0.116 (0.03)"

write.table(yearz2_full_table, file = "year_by_zone.csv", row.names = TRUE,
            col.names = NA, sep = ",")

Without Cape Sebastian or Pistol River South Fork

noseb <- popsub(setpop(ramdat, ~ZONE2), blacklist = "HunterCr")

system.time(noseb_year_z2   <- poppr.amova(noseb, ~Year/ZONE2))

## 
##  No missing values detected.

##    user  system elapsed 
##   0.776   0.000   0.778

system.time(noseb_year_z2cc <- poppr.amova(noseb, ~Year/ZONE2, clonecorrect = TRUE))

## 
##  No missing values detected.

##    user  system elapsed 
##   0.892   0.000   0.892

set.seed(9001)
system.time(noseb_year_z2.t <- randtest(noseb_year_z2, nrepet = 999))

##    user  system elapsed 
## 114.244   0.664 114.935

plot(noseb_year_z2.t)

set.seed(9001)
system.time(noseb_year_z2cc.t <- randtest(noseb_year_z2cc, nrepet = 999))

##    user  system elapsed 
##  15.653   0.000  15.658

plot(noseb_year_z2cc.t)

date()

## [1] "Fri Oct 24 16:23:52 2014"

noseb <- popsub(setpop(ramdat, ~ZONE2), blacklist = "HunterCr")

system.time(noseb_z2_year   <- poppr.amova(noseb, ~ZONE2/Year))

## 
##  No missing values detected.

##    user  system elapsed 
##   0.768   0.052   0.821

system.time(noseb_z2_yearcc <- poppr.amova(noseb, ~ZONE2/Year, clonecorrect = TRUE))

## 
##  No missing values detected.

##    user  system elapsed 
##   0.900   0.016   0.915

set.seed(9001)
system.time(noseb_z2_year.t <- randtest(noseb_z2_year, nrepet = 999))

##    user  system elapsed 
## 113.627   1.380 115.031

plot(noseb_z2_year.t)

set.seed(9001)
system.time(noseb_z2_yearcc.t <- randtest(noseb_z2_yearcc, nrepet = 999))

##    user  system elapsed 
##  15.425   0.200  15.629

plot(noseb_z2_yearcc.t)

date()

## [1] "Fri Oct 24 16:45:05 2014"

nosebpr <- selPopSize(noseb, n = 10)

system.time(nosebpr_year_z2   <- poppr.amova(nosebpr, ~Year/ZONE2))

## 
##  No missing values detected.

##    user  system elapsed 
##   0.764   0.004   0.767

system.time(nosebpr_year_z2cc <- poppr.amova(nosebpr, ~Year/ZONE2, clonecorrect = TRUE))

## 
##  No missing values detected.

##    user  system elapsed 
##   0.852   0.024   0.875

set.seed(9001)
system.time(nosebpr_year_z2.t <- randtest(nosebpr_year_z2, nrepet = 999))

##    user  system elapsed 
## 112.171   1.472 113.661

plot(nosebpr_year_z2.t)

set.seed(9001)
system.time(nosebpr_year_z2cc.t <- randtest(nosebpr_year_z2cc, nrepet = 999))

##    user  system elapsed 
##  15.269   0.216  15.489

plot(nosebpr_year_z2cc.t)

nosebpr <- selPopSize(noseb, n = 10)

system.time(nosebpr_z2_year   <- poppr.amova(nosebpr, ~ZONE2/Year))

## 
##  No missing values detected.

##    user  system elapsed 
##   0.780   0.012   0.793

system.time(nosebpr_z2_yearcc <- poppr.amova(nosebpr, ~ZONE2/Year, clonecorrect = TRUE))

## 
##  No missing values detected.

##    user  system elapsed 
##   0.860   0.008   0.867

set.seed(9001)
system.time(nosebpr_z2_year.t <- randtest(nosebpr_z2_year, nrepet = 999))

##    user  system elapsed 
## 110.799   1.376 112.190

plot(nosebpr_z2_year.t)

set.seed(9001)
system.time(nosebpr_z2_yearcc.t <- randtest(nosebpr_z2_yearcc, nrepet = 999))

##    user  system elapsed 
##  15.165   0.160  15.326

plot(nosebpr_z2_yearcc.t)

Excluding individual variance and Cape Sebastian/PRSF

Year/Region

nosebpr <- resetMLG(nosebpr)

(nosebpr_WF_year_z2   <- poppr.amova(nosebpr, ~Year/ZONE2, within = FALSE))

## 
##  No missing values detected.

## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
## 
## $results
##                              Df Sum Sq Mean Sq
## Between Year                 10  48.86  4.8862
## Between samples Within Year  10  54.55  5.4547
## Within samples              422 272.72  0.6463
## Total                       442 376.13  0.8510
## 
## $componentsofcovariance
##                                           Sigma      %
## Variations  Between Year                -0.1154 -13.31
## Variations  Between samples Within Year  0.3359  38.75
## Variations  Within samples               0.6463  74.56
## Total variations                         0.8668 100.00
## 
## $statphi
##                       Phi
## Phi-samples-total  0.2544
## Phi-samples-Year   0.3420
## Phi-Year-total    -0.1331

(nosebpr_WF_year_z2cc <- poppr.amova(nosebpr, ~Year/ZONE2, clonecorrect = TRUE, 
                          within = FALSE))

## 
##  No missing values detected.

## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
## 
## $results
##                              Df Sum Sq Mean Sq
## Between Year                 10  18.95   1.895
## Between samples Within Year  10  16.47   1.647
## Within samples              125 138.04   1.104
## Total                       145 173.45   1.196
## 
## $componentsofcovariance
##                                           Sigma       %
## Variations  Between Year                0.02078   1.718
## Variations  Between samples Within Year 0.08416   6.960
## Variations  Within samples              1.10429  91.322
## Total variations                        1.20923 100.000
## 
## $statphi
##                       Phi
## Phi-samples-total 0.08678
## Phi-samples-Year  0.07082
## Phi-Year-total    0.01718

set.seed(9001)
system.time(nosebpr_WF_year_z2.t <- randtest(nosebpr_WF_year_z2, nrepet = 9999))

##    user  system elapsed 
##  27.845   0.185  28.338

plot(nosebpr_WF_year_z2.t)

plot of chunk nosebpr_AMOVA_WF

set.seed(9001)
system.time(nosebpr_WF_year_z2cc.t <- randtest(nosebpr_WF_year_z2cc, nrepet = 9999))

##    user  system elapsed 
##  16.948   0.205  19.099

plot(nosebpr_WF_year_z2cc.t)

plot of chunk nosebpr_AMOVA_WF

date()

## [1] "Fri Nov 14 10:34:00 2014"

Region/Year

namehierarchy(nosebpr) <- ~Year/ZONE2/ZONE1

(nosebpr_WF_z2_year   <- poppr.amova(nosebpr, ~ZONE2/Year, within = FALSE))

## 
##  No missing values detected.

## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
## 
## $results
##                               Df Sum Sq Mean Sq
## Between ZONE2                  4  80.91 20.2284
## Between samples Within ZONE2  16  22.49  1.4059
## Within samples               422 272.72  0.6463
## Total                        442 376.13  0.8510
## 
## $componentsofcovariance
##                                            Sigma       %
## Variations  Between ZONE2                0.26437  27.599
## Variations  Between samples Within ZONE2 0.04727   4.935
## Variations  Within samples               0.64627  67.467
## Total variations                         0.95791 100.000
## 
## $statphi
##                       Phi
## Phi-samples-total 0.32533
## Phi-samples-ZONE2 0.06816
## Phi-ZONE2-total   0.27599

(nosebpr_WF_z2_yearcc <- poppr.amova(nosebpr, ~ZONE2/Year, clonecorrect = TRUE, 
                          within = FALSE))

## 
##  No missing values detected.

## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
## 
## $results
##                               Df Sum Sq Mean Sq
## Between ZONE2                  4  18.21   4.553
## Between samples Within ZONE2  16  17.20   1.075
## Within samples               125 138.04   1.104
## Total                        145 173.45   1.196
## 
## $componentsofcovariance
##                                              Sigma        %
## Variations  Between ZONE2                 0.138799  11.2119
## Variations  Between samples Within ZONE2 -0.005119  -0.4135
## Variations  Within samples                1.104290  89.2016
## Total variations                          1.237970 100.0000
## 
## $statphi
##                         Phi
## Phi-samples-total  0.107984
## Phi-samples-ZONE2 -0.004657
## Phi-ZONE2-total    0.112119

set.seed(9001)
system.time(nosebpr_WF_z2_year.t <- randtest(nosebpr_WF_z2_year, nrepet = 9999))

##    user  system elapsed 
##  26.774   0.203  27.382

plot(nosebpr_WF_z2_year.t)

plot of chunk nosebpr_AMOVAzy_WF

set.seed(9001)
system.time(nosebpr_WF_z2_yearcc.t <- randtest(nosebpr_WF_z2_yearcc, nrepet = 9999))

##    user  system elapsed 
##  14.217   0.091  14.394

plot(nosebpr_WF_z2_yearcc.t)

plot of chunk nosebpr_AMOVAzy_WF

date()

## [1] "Fri Nov 14 10:34:43 2014"

nosebpr_z2year_table   <- make_amova_table(nosebpr_WF_z2_year, nosebpr_WF_z2_year.t)
nosebpr_z2yearcc_table <- make_amova_table(nosebpr_WF_z2_yearcc, nosebpr_WF_z2_yearcc.t)
(nosebpr_z2year_full_table <- make_amova_printable(nosebpr_z2year_table, nosebpr_z2yearcc_table))

##                               statistic
## levels                         d.f.        Sum of Squares
##   Between ZONE2                "4 (4)"     "80.9 (18.2)" 
##   Between samples Within ZONE2 "16 (16)"   "22.5 (17.2)" 
##   Within samples               "422 (125)" "273 (138)"   
##                               statistic
## levels                         Percent variation P              
##   Between ZONE2                "27.6 (11.2)"     "1e-04 (1e-04)"
##   Between samples Within ZONE2 "4.93 (-0.413)"   "1e-04 (0.445)"
##   Within samples               "67.5 (89.2)"     "1e-04 (1e-04)"
##                               statistic
## levels                         Phi statistic      
##   Between ZONE2                "0.325 (0.108)"    
##   Between samples Within ZONE2 "0.0682 (-0.00466)"
##   Within samples               "0.276 (0.112)"

write.table(nosebpr_z2year_full_table, file = "nosebpr_zone_by_year.csv", row.names = TRUE,
            col.names = NA, sep = ",")

nosebpr_yearz2_table   <- make_amova_table(nosebpr_WF_year_z2, nosebpr_WF_year_z2.t)
nosebpr_yearz2cc_table <- make_amova_table(nosebpr_WF_year_z2cc, nosebpr_WF_year_z2cc.t)
(nosebpr_yearz2_full_table <- make_amova_printable(nosebpr_yearz2_table, nosebpr_yearz2cc_table))

##                              statistic
## levels                        d.f.        Sum of Squares Percent variation
##   Between Year                "10 (10)"   "48.9 (18.9)"  "-13.3 (1.72)"   
##   Between samples Within Year "10 (10)"   "54.5 (16.5)"  "38.8 (6.96)"    
##   Within samples              "422 (125)" "273 (138)"    "74.6 (91.3)"    
##                              statistic
## levels                        P               Phi statistic    
##   Between Year                "0.796 (0.293)" "0.254 (0.0868)" 
##   Between samples Within Year "1e-04 (4e-04)" "0.342 (0.0708)" 
##   Within samples              "1e-04 (1e-04)" "-0.133 (0.0172)"

write.table(nosebpr_yearz2_full_table, file = "nosebpr_year_by_zone.csv", row.names = TRUE,
            col.names = NA, sep = ",")

Without the coast, cape sebastian, or PRSF

Year/Region

nosebprc <- popsub(setpop(nosebpr, ~ZONE2), blacklist = "Coast")
nosebprc <- resetMLG(nosebprc)

(nosebprc_WF_year_z2   <- poppr.amova(nosebprc, ~Year/ZONE2, within = FALSE))

## 
##  No missing values detected.

## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
## 
## $results
##                              Df Sum Sq Mean Sq
## Between Year                  7  41.20  5.8852
## Between samples Within Year   7  37.96  5.4226
## Within samples              394 253.02  0.6422
## Total                       408 332.17  0.8141
## 
## $componentsofcovariance
##                                           Sigma      %
## Variations  Between Year                -0.1232 -14.83
## Variations  Between samples Within Year  0.3116  37.52
## Variations  Within samples               0.6422  77.31
## Total variations                         0.8306 100.00
## 
## $statphi
##                       Phi
## Phi-samples-total  0.2269
## Phi-samples-Year   0.3267
## Phi-Year-total    -0.1483

(nosebprc_WF_year_z2cc <- poppr.amova(nosebprc, ~Year/ZONE2, clonecorrect = TRUE, 
                          within = FALSE))

## 
##  No missing values detected.

## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
## 
## $results
##                              Df Sum Sq Mean Sq
## Between Year                  7  14.30   2.043
## Between samples Within Year   7  10.03   1.433
## Within samples              112 125.39   1.120
## Total                       126 149.72   1.188
## 
## $componentsofcovariance
##                                           Sigma       %
## Variations  Between Year                0.03813   3.168
## Variations  Between samples Within Year 0.04603   3.824
## Variations  Within samples              1.11958  93.008
## Total variations                        1.20375 100.000
## 
## $statphi
##                       Phi
## Phi-samples-total 0.06992
## Phi-samples-Year  0.03949
## Phi-Year-total    0.03168

set.seed(9001)
system.time(nosebprc_WF_year_z2.t <- randtest(nosebprc_WF_year_z2, nrepet = 9999))

##    user  system elapsed 
##  18.395   0.074  18.592

plot(nosebprc_WF_year_z2.t)

plot of chunk nosebprc_AMOVA_WF

set.seed(9001)
system.time(nosebprc_WF_year_z2cc.t <- randtest(nosebprc_WF_year_z2cc, nrepet = 9999))

##    user  system elapsed 
##   9.906   0.040   9.966

plot(nosebprc_WF_year_z2cc.t)

plot of chunk nosebprc_AMOVA_WF

date()

## [1] "Sun Nov 16 09:51:22 2014"

Region/Year

namehierarchy(nosebprc) <- ~Year/ZONE2/ZONE1

(nosebprc_WF_z2_year   <- poppr.amova(nosebprc, ~ZONE2/Year, within = FALSE))

## 
##  No missing values detected.

## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
## 
## $results
##                               Df Sum Sq Mean Sq
## Between ZONE2                  3  62.36 20.7878
## Between samples Within ZONE2  11  16.79  1.5265
## Within samples               394 253.02  0.6422
## Total                        408 332.17  0.8141
## 
## $componentsofcovariance
##                                            Sigma       %
## Variations  Between ZONE2                0.24461  26.354
## Variations  Between samples Within ZONE2 0.04136   4.456
## Variations  Within samples               0.64217  69.189
## Total variations                         0.92814 100.000
## 
## $statphi
##                       Phi
## Phi-samples-total 0.30811
## Phi-samples-ZONE2 0.06051
## Phi-ZONE2-total   0.26354

(nosebprc_WF_z2_yearcc <- poppr.amova(nosebprc, ~ZONE2/Year, clonecorrect = TRUE, 
                          within = FALSE))

## 
##  No missing values detected.

## $call
## ade4::amova(samples = xtab, distances = xdist, structures = xstruct)
## 
## $results
##                               Df Sum Sq Mean Sq
## Between ZONE2                  3  11.59   3.863
## Between samples Within ZONE2  11  12.74   1.158
## Within samples               112 125.39   1.120
## Total                        126 149.72   1.188
## 
## $componentsofcovariance
##                                             Sigma       %
## Variations  Between ZONE2                0.103175   8.400
## Variations  Between samples Within ZONE2 0.005491   0.447
## Variations  Within samples               1.119584  91.153
## Total variations                         1.228250 100.000
## 
## $statphi
##                       Phi
## Phi-samples-total 0.08847
## Phi-samples-ZONE2 0.00488
## Phi-ZONE2-total   0.08400

set.seed(9001)
system.time(nosebprc_WF_z2_year.t <- randtest(nosebprc_WF_z2_year, nrepet = 9999))

##    user  system elapsed 
##  17.590   0.057  17.679

plot(nosebprc_WF_z2_year.t)

plot of chunk nosebprc_AMOVAzy_WF

set.seed(9001)
system.time(nosebprc_WF_z2_yearcc.t <- randtest(nosebprc_WF_z2_yearcc, nrepet = 9999))

##    user  system elapsed 
##   9.205   0.036   9.258

plot(nosebprc_WF_z2_yearcc.t)

plot of chunk nosebprc_AMOVAzy_WF

date()

## [1] "Sun Nov 16 09:51:53 2014"

nosebprc_z2year_table   <- make_amova_table(nosebprc_WF_z2_year, nosebprc_WF_z2_year.t)
nosebprc_z2yearcc_table <- make_amova_table(nosebprc_WF_z2_yearcc, nosebprc_WF_z2_yearcc.t)
(nosebprc_z2year_full_table <- make_amova_printable(nosebprc_z2year_table, nosebprc_z2yearcc_table))

##                               statistic
## levels                         d.f.        Sum of Squares
##   Between ZONE2                "3 (3)"     "62.4 (11.6)" 
##   Between samples Within ZONE2 "11 (11)"   "16.8 (12.7)" 
##   Within samples               "394 (112)" "253 (125)"   
##                               statistic
## levels                         Percent variation P               
##   Between ZONE2                "26.4 (8.4)"      "2e-04 (0.0012)"
##   Between samples Within ZONE2 "4.46 (0.447)"    "1e-04 (0.362)" 
##   Within samples               "69.2 (91.2)"     "1e-04 (1e-04)" 
##                               statistic
## levels                         Phi statistic     
##   Between ZONE2                "0.308 (0.0885)"  
##   Between samples Within ZONE2 "0.0605 (0.00488)"
##   Within samples               "0.264 (0.084)"

write.table(nosebprc_z2year_full_table, file = "nosebprc_zone_by_year.csv", row.names = TRUE,
            col.names = NA, sep = ",")

nosebprc_yearz2_table   <- make_amova_table(nosebprc_WF_year_z2, nosebprc_WF_year_z2.t)
nosebprc_yearz2cc_table <- make_amova_table(nosebprc_WF_year_z2cc, nosebprc_WF_year_z2cc.t)
(nosebprc_yearz2_full_table <- make_amova_printable(nosebprc_yearz2_table, nosebprc_yearz2cc_table))

##                              statistic
## levels                        d.f.        Sum of Squares Percent variation
##   Between Year                "7 (7)"     "41.2 (14.3)"  "-14.8 (3.17)"   
##   Between samples Within Year "7 (7)"     "38 (10)"      "37.5 (3.82)"    
##   Within samples              "394 (112)" "253 (125)"    "77.3 (93)"      
##                              statistic
## levels                        P                Phi statistic    
##   Between Year                "0.843 (0.244)"  "0.227 (0.0699)" 
##   Between samples Within Year "1e-04 (0.0265)" "0.327 (0.0395)" 
##   Within samples              "1e-04 (3e-04)"  "-0.148 (0.0317)"

write.table(nosebprc_yearz2_full_table, file = "nosebprc_year_by_zone.csv", row.names = TRUE,
            col.names = NA, sep = ",")

Pairwise AMOVA Between Regions

While presenting the AMOVA results, Ebba brought up a very good point that the more recently introduced regions might be biasing the AMOVA results. To assess which regions might have differentiation in a pairwise fashion, a pairwise AMOVA is being run on regions comparing both regions separately per year.

Running the AMOVA by zone

zonePairs    <- pairwise_amova(setpop(ramdat, ~ZONE2), ~ZONE2/Year, quiet = TRUE)
zonePairs.cc <- pairwise_amova(setpop(ramdat, ~ZONE2), ~ZONE2/Year, 
                                quiet = TRUE, clonecorrect = TRUE)

set.seed(9001) #
system.time(zonePairs.t    <- pairwise_amova_test(zonePairs, nrepet = 999))

## JHallCr : NFChetHigh

## JHallCr : Coast

## JHallCr : HunterCr

## JHallCr : Winchuck

## JHallCr : ChetcoMain

## JHallCr : PistolRSF

## NFChetHigh : Coast

## NFChetHigh : HunterCr

## NFChetHigh : Winchuck

## NFChetHigh : ChetcoMain

## NFChetHigh : PistolRSF

## Coast : HunterCr

## Coast : Winchuck

## Coast : ChetcoMain

## Coast : PistolRSF

## HunterCr : Winchuck

## HunterCr : ChetcoMain

## HunterCr : PistolRSF

## Winchuck : ChetcoMain

## Winchuck : PistolRSF

## ChetcoMain : PistolRSF

##    user  system elapsed 
## 297.830   0.140 298.174

set.seed(9001)
system.time(zonePairs.cc.t <- pairwise_amova_test(zonePairs.cc, nrepet = 999))

## JHallCr : NFChetHigh

## JHallCr : Coast

## JHallCr : HunterCr

## JHallCr : Winchuck

## JHallCr : ChetcoMain

## JHallCr : PistolRSF

## NFChetHigh : Coast

## NFChetHigh : HunterCr

## NFChetHigh : Winchuck

## NFChetHigh : ChetcoMain

## NFChetHigh : PistolRSF

## Coast : HunterCr

## Coast : Winchuck

## Coast : ChetcoMain

## Coast : PistolRSF

## HunterCr : Winchuck

## HunterCr : ChetcoMain

## HunterCr : PistolRSF

## Winchuck : ChetcoMain

## Winchuck : PistolRSF

## ChetcoMain : PistolRSF

##    user  system elapsed 
##  42.071   0.144  42.224

Running the AMOVA by year

yearPairs    <- pairwise_amova(setpop(ramdat, ~Year), ~Year/ZONE2, quiet = TRUE)
yearPairs.cc <- pairwise_amova(setpop(ramdat, ~Year), ~Year/ZONE2, 
                                quiet = TRUE, clonecorrect = TRUE)

set.seed(9001)
system.time(yearPairs.t    <- pairwise_amova_test(yearPairs, nrepet = 999))

## 2001 : 2002

## 2001 : 2003

## 2001 : 2004

## 2001 : 2011

## 2001 : 2012

## 2001 : 2013

## 2001 : 2014

## 2002 : 2003

## 2002 : 2004

## 2002 : 2011

## 2002 : 2012

## 2002 : 2013

## 2002 : 2014

## 2003 : 2004

## 2003 : 2011

## 2003 : 2012

## 2003 : 2013

## 2003 : 2014

## 2004 : 2011

## 2004 : 2012

## 2004 : 2013

## 2004 : 2014

## 2011 : 2012

## 2011 : 2013

## 2011 : 2014

## 2012 : 2013

## 2012 : 2014

## 2013 : 2014

##    user  system elapsed 
## 276.694   1.860 279.005

set.seed(9001)
system.time(yearPairs.cc.t <- pairwise_amova_test(yearPairs.cc, nrepet = 999))

## 2001 : 2002

## 2001 : 2003

## 2001 : 2004

## 2001 : 2011

## 2001 : 2012

## 2001 : 2013

## 2001 : 2014

## 2002 : 2003

## 2002 : 2004

## 2002 : 2011

## 2002 : 2012

## 2002 : 2013

## 2002 : 2014

## 2003 : 2004

## 2003 : 2011

## 2003 : 2012

## 2003 : 2013

## 2003 : 2014

## 2004 : 2011

## 2004 : 2012

## 2004 : 2013

## 2004 : 2014

## 2011 : 2012

## 2011 : 2013

## 2011 : 2014

## 2012 : 2013

## 2012 : 2014

## 2013 : 2014

##    user  system elapsed 
##  44.011   0.100  44.107

Pairwise Heatmaps

By region

names(zonePairs.t)    <- names(zonePairs)
names(zonePairs.cc.t) <- names(zonePairs.cc)
parray <- array(dim = c(length(zonePairs.t[[1]]$pvalue),
                        length(zonePairs.t),
                        2),
                dimnames = list(Test = zonePairs.t[[1]]$names,
                                Pair = names(zonePairs.t),
                                data = c("full", "clone-censored"))
                )
parray[, , "full"]           <- vapply(zonePairs.t, getpval,
                                       zonePairs.t[[1]]$pvalue)

parray[, , "clone-censored"] <- vapply(zonePairs.cc.t, getpval,
                                       zonePairs.cc.t[[1]]$pvalue)

phiarray <- array(dim = c(length(getphi(zonePairs[[1]])),
                        length(zonePairs),
                        2),
                dimnames = list(Test = rownames(zonePairs[[1]]$statphi),
                                Pair = names(zonePairs),
                                data = c("full", "clone-censored"))
                )
phiarray[, , "full"]           <- vapply(zonePairs, getphi,
                                         getphi(zonePairs[[1]]))

phiarray[, , "clone-censored"] <- vapply(zonePairs.cc, getphi,
                                         getphi(zonePairs.cc[[1]]))



p.df <- tbl_df(melt(parray[3:4, , ])) %>% 
  separate(Pair, c("PairOne", "PairTwo"), sep = " : ")

# Fixing the data so it appaears as an upper-triangle
setpop(ramdat) <- ~ZONE2
p.df$PairOne <- refactor(p.df$PairOne, pop(ramdat))
p.df$PairTwo <- refactor(p.df$PairTwo, pop(ramdat))
# p.df$PairOne <- unique(pop(ramdat))[match(p.df$PairOne, unique(pop(ramdat)))]
# p.df$PairTwo <- unique(pop(ramdat))[match(p.df$PairTwo, unique(pop(ramdat)))]
levels(p.df$Test) <- paste("variations between", c("year within region", "regions"))

# Plotting
ggplot(p.df, aes(x = PairTwo, y = PairOne, fill = -log(value))) + 
  geom_tile() + 
  geom_text(aes(label = value, color = value)) +
  facet_grid(data ~ Test) + 
  scale_fill_continuous(low = "black", high = "orange", name = "p value\n",
                        breaks = c(0.693, 3, 4.6, 6.7),
                        labels = round(exp(-c(0.693, 3, 4.6, 6.8)), 3)
                        ) +
  scale_color_continuous(low = "black", high = "orange", guide = "none") +
  theme_classic() + myTheme + labs(list(x = NULL, y = NULL))

## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_text).
## Warning: Removed 1 rows containing missing values (geom_text).

plot of chunk plot_amova_pvals

# Creating a matrix to store Phi values
phi <- matrix(nrow = nlevels(p.df$PairOne), ncol = nlevels(p.df$PairOne), 
              dimnames = list(levels(p.df$PairOne), levels(p.df$PairOne)))
p   <- phi
phi[lower.tri(phi)] <- phiarray[4, , 2] # Clone corrected (upper triangle)
phi <- t(phi) # Transposition to final form
phi[lower.tri(phi)] <- phiarray[4, , 1] # Raw data (lower triangle)
phi

##             JHallCr NFChetHigh    Coast HunterCr Winchuck ChetcoMain
## JHallCr          NA    0.01399 0.041699  0.08363  0.03166   0.017933
## NFChetHigh  0.02571         NA 0.024833  0.07175  0.02982   0.018233
## Coast       0.06324    0.04996       NA  0.04240  0.05376   0.002047
## HunterCr    0.18133    0.14962 0.082444       NA  0.11903   0.019447
## Winchuck    0.06174    0.06309 0.071333  0.18991       NA   0.026352
## ChetcoMain  0.08349    0.08168 0.002977 -0.02980  0.07668         NA
## PistolRSF  -0.01156    0.01162 0.035642      NaN  0.03453  -0.061292
##            PistolRSF
## JHallCr    -0.031101
## NFChetHigh -0.024162
## Coast       0.014868
## HunterCr         NaN
## Winchuck   -0.003834
## ChetcoMain -0.037686
## PistolRSF         NA

p[lower.tri(p)]     <- parray[4, , 2]
p <- t(p)
p[lower.tri(p)]     <- parray[4, , 1]
p

##            JHallCr NFChetHigh Coast HunterCr Winchuck ChetcoMain PistolRSF
## JHallCr         NA      0.963 1.000    1.000    1.000      0.977     0.001
## NFChetHigh   0.996         NA 0.986    1.000    0.934      1.000     0.205
## Coast        1.000      0.977    NA    0.866    1.000      0.680     0.730
## HunterCr     1.000      1.000 1.000       NA    1.000      1.000        NA
## Winchuck     1.000      0.938 1.000    1.000       NA      1.000     0.667
## ChetcoMain   1.000      0.866 0.710    0.334    1.000         NA     0.341
## PistolRSF    0.001      0.809 0.720       NA    1.000      0.353        NA

write.table(phi, file = "phi_by_zone.csv", sep = ",", row.names = TRUE,
            col.names = NA)

By year

names(yearPairs.t)    <- names(yearPairs)
names(yearPairs.cc.t) <- names(yearPairs.cc)
parray_year <- array(dim = c(length(yearPairs.t[[1]]$pvalue),
                        length(yearPairs.t),
                        2),
                dimnames = list(Test = yearPairs.t[[1]]$names,
                                Pair = names(yearPairs.t),
                                data = c("full", "clone-censored"))
                )

parray_year[, , "full"]           <- vapply(yearPairs.t, getpval, 
                                            yearPairs.t[[1]]$pvalue)

parray_year[, , "clone-censored"] <- vapply(yearPairs.cc.t, getpval,
                                            yearPairs.cc.t[[1]]$pvalue)

phiarray_year <- array(dim = c(length(getphi(yearPairs[[1]])),
                        length(yearPairs),
                        2),
                dimnames = list(Test = rownames(yearPairs[[1]]$statphi),
                                Pair = names(yearPairs),
                                data = c("full", "clone-censored"))
                )
phiarray_year[, , "full"]           <- vapply(yearPairs, getphi,
                                         getphi(yearPairs[[1]]))

phiarray_year[, , "clone-censored"] <- vapply(yearPairs.cc, getphi,
                                         getphi(yearPairs.cc[[1]]))

py.df <- tbl_df(melt(parray_year[3:4, , ])) %>% 
  separate(Pair, c("PairOne", "PairTwo"), sep = " : ")

# Fixing the data so it appaears as an upper-triangle
setpop(ramdat) <- ~Year
py.df$PairOne <- refactor(py.df$PairOne, pop(ramdat))
py.df$PairTwo <- refactor(py.df$PairTwo, pop(ramdat))
# py.df$PairOne <- unique(pop(ramdat))[match(py.df$PairOne, unique(pop(ramdat)))]
# py.df$PairTwo <- unique(pop(ramdat))[match(py.df$PairTwo, unique(pop(ramdat)))]
levels(py.df$Test) <- paste("variations between", c("regions within years", "years"))

# Plotting
ggplot(py.df, aes(x = PairTwo, y = PairOne, fill = -log(value))) + 
  geom_tile() + 
  geom_text(aes(label = value, color = value)) +
  facet_grid(data ~ Test) + 
  scale_fill_continuous(low = "black", high = "orange", name = "p value\n",
                        breaks = c(0.693, 3, 4.6, 6.7),
                        labels = round(exp(-c(0.693, 3, 4.6, 6.8)), 3)
                        ) +
  scale_color_continuous(low = "black", high = "orange", guide = "none") +
  theme_classic() + myTheme + labs(list(x = NULL, y = NULL))

## Warning: Removed 30 rows containing missing values (geom_text).
## Warning: Removed 30 rows containing missing values (geom_text).
## Warning: Removed 30 rows containing missing values (geom_text).
## Warning: Removed 30 rows containing missing values (geom_text).

plot of chunk plot_amova_pvals_year

phi_year <- matrix(nrow = nlevels(py.df$PairOne), ncol = nlevels(py.df$PairOne), 
              dimnames = list(levels(py.df$PairOne), levels(py.df$PairOne)))

phi_year[lower.tri(phi_year)] <- phiarray_year[4, , 2]
phi_year <- t(phi_year)
phi_year[lower.tri(phi_year)] <- phiarray_year[4, , 1]
phi_year

##           2001     2002      2003     2004 2005 2006 2010     2011
## 2001        NA      NaN  0.041386      NaN    0    0    0 -0.02217
## 2002       NaN       NA  0.033406      NaN    0    0    0 -0.01252
## 2003 -0.018524 -0.01675        NA 0.043358    0    0    0  0.02827
## 2004       NaN      NaN -0.004773       NA    0    0    0  0.01991
## 2005  0.000000  0.00000  0.000000 0.000000   NA    0    0  0.00000
## 2006  0.000000  0.00000  0.000000 0.000000    0   NA    0  0.00000
## 2010  0.000000  0.00000  0.000000 0.000000    0    0   NA  0.00000
## 2011  0.031110  0.02759  0.075251 0.052300    0    0    0       NA
## 2012 -0.016712 -0.01932 -0.002365 0.018806    0    0    0  0.08346
## 2013 -0.013666 -0.02135 -0.015477 0.003758    0    0    0  0.09174
## 2014 -0.005429 -0.01381 -0.004931 0.011101    0    0    0  0.06162
##            2012       2013       2014
## 2001 -0.0182409  1.669e-03  0.0059683
## 2002 -0.0003216 -2.695e-05 -0.0042011
## 2003  0.0175266  1.206e-02  0.0146585
## 2004  0.0405752  4.187e-02  0.0311115
## 2005  0.0000000  0.000e+00  0.0000000
## 2006  0.0000000  0.000e+00  0.0000000
## 2010  0.0000000  0.000e+00  0.0000000
## 2011  0.0161974  1.903e-02 -0.0006563
## 2012         NA -1.181e-02 -0.0059601
## 2013 -0.0041298         NA -0.0034361
## 2014  0.0074928 -1.831e-02         NA

write.table(phi_year, file = "phi_by_year.csv", sep = ",", row.names = TRUE,
            col.names = NA)

Diversity_statistics

Zhian N. Kamvar

October 15, 2014

Purpose

Required packages and data

Functions

Genotypic Diversity

Confidence intervals for genotypic diversity

Allelic and Genotypic diversity per locus

Tables

Mean Allelic Diversity by Region

Mean Allelic Diversity by Year

Total Allelic Diversity

AMOVA

Including within individual variance

Year/Region

Region/Year

Excluding within individual variance

Year/Region

Region/Year

Writing tables

Without Cape Sebastian or Pistol River South Fork

Excluding individual variance and Cape Sebastian/PRSF

Year/Region

Region/Year

Without the coast, cape sebastian, or PRSF

Year/Region

Region/Year

Pairwise AMOVA Between Regions

Running the AMOVA by zone

Running the AMOVA by year

Pairwise Heatmaps

By region

By year